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This paper introduces a Monte Carlo method for maximum hke- 
hhood inference in the context of discretely observed diffusion pro- 
cesses. The method gives unbiased and a.s. continuous estimators of 
the likelihood function for a family of diffusion models and its per- 
formance in numerical examples is computationally efficient. It uses 
a recently developed technique for the exact simulation of diffusions, 
and involves no discretization error. We show that, under regularity 
conditions, the Monte Carlo MLE converges a.s. to the true MLE. 
For datasize n — > oo, we show that the number of Monte Carlo itera- 
tions should be tuned as 0{n^^^) and we demonstrate the consistency 
properties of the Monte Carlo MLE as an estimator of the true pa- 
rameter value. 

1. Introduction. We introduce a Monte Carlo method for maximum like- 
lihood inference in the context of discretely observed diffusion processes. The 
method gives unbiased and a.s. continuous estimates of the likelihood func- 
tion, which converge uniformly in the parameters to the likelihood function 
as the Monte Carlo sample size N increases. Additionally, for increasing 
datasize n, the asymptotically optimal algorithm corresponds to selecting 
iV = 0(ni/2). 

Consider scalar time-homogeneous diffusion processes, defined by stochas- 
tic differential equations (SDEs) of the type 

(1) dVs = b{Vs;e)ds + a{Vs;e)dBs, € V C R, 

where B is a Brownian motion. The drift and the diffusion coefficient, b{-;6) 
and cr(-; 0) respectively, are assumed to be known up to a vector of parame- 
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ters d & Q C R*^. Multivariate extensions of this work are addressed in Sec- 
tion 7. The process is assumed to be observed without error at a cohection 
of time instances, 

{Vt„Vt„...,VtJ, = 1 < i 1 < • • • < i n , 

and the statistical challenge is to explore characteristics, such as the maxi- 
mum and the level sets, of the likelihood function 

n 

1=1 

where Ati = ti — ti-i and 

(2) pt{v,w;0) :=F[Vtedw\Vo = v;9]/dw, t>0,w,veV, 

is the transition density of (1). It is well documented (see, e.g., [37] for a 
recent review) that inference about diffusion models is complicated by the 
unavailability of the transition density (2), except for limited cases. Thus, in- 
ference is carried out using either nonlikelihood approaches, such as estimat- 
ing equations [10], efficient method of moments [20] and indirect inference 
[25], or approximate likelihood-based approaches, such as computationally 
intensive Markov chain Monte Carlo imputation methods [17, 18, 35] and 
methods which approximate analytically the transition density [1]. Also, 
approximate Monte Carlo maximum likelihood approaches have been sug- 
gested, most notably by [32] and [16]. 

The method described in this paper, termed the Simultaneous Acceptance 
Method (SAM), estimates each likelihood contribution Li{-) independently. 
To simplify the presentation, let L(-) denote the likelihood contribution of 
an arbitrary pair of consecutive data points Vs = v, Vg+t = w, that is, L{9) = 
pt{v,w;9). SAM generates a random function L{'E,6), 9 £&, where H is a 
random element independent of 9, such that for any fixed 9 gQ, E[L(H, 9)] = 
L{9), and w.p.l, 9 L(E,9) is continuous. We estimate L(-) by the Monte 
Carlo functional averages 

(3) L'^ {■) := — ''^L{=.\-), =.^,...,=. independent copies of r.. 

Having obtained estimates Lf{-) of Li{-) using (3) independently for each 
i = 1, . . . ,n, we estimate the likelihood function £«(•) by the product 

n 

£;r(-):=n^f(-)- 

i=l 

We demonstrate that our Monte Carlo estimator £^ (•) is computation- 
ally efficient and we detail its theoretical properties. From the computational 
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perspective, the random element consists of standard exponential and Gaus- 
sian variables and can be easily simulated, and L{E,9) is of a simple calcu- 
lable form. Moreover, the a.s. continuity of (•) facilitates the efficient 
implementation of optimization routines for locating its maximizer. From 
the theoretical perspective, £^ (9) is an unbiased estimator of -Cn(^) with 
finite moments for any 9 £ @. More importantly, due to the a.s. continuity 
of £^ (•), we can resort to the Strong Law of Large Numbers (SLLN) in 
Banach spaces, to establish that w.p.l £^ (6) converges to Sln{9) uniformly 
in 9, as — > oo. Uniform convergence to the likelihood implies, under mild 
conditions, convergence of the maximizers to the MLE, say 9n- 

(4) w.p.l, 9^ ■.= aigmax£,^ {9)^ 9n asN^oo. 

The rate of convergence is shown to be 0{N^^'^). Also, additional statistics, 
for example, profile likelihoods and level sets, can be derived as appropriate 
limits of the corresponding characteristics of the Monte Carlo estimate of 
the likelihood. 

We also investigate the properties of the Monte Carlo method for increas- 
ing datasize n. The Monte Carlo MLE 9^ is a consistent estimator of the true 
parameter value, say, when N — > oo, n — > oo. The optimal algorithm, in 
terms of computational cost as n — > oo, is obtained when N = Nn = 0{n^^'^), 
whence 7i^/^{9^ -9o) converges in law to a Gaussian random variable. 

The construction of the random function L{E,-) is based on a recently 
developed retrospective rejection sampling algorithm called the Exact Algo- 
rithm (EA). EA was introduced in [9] and [7]; it returns a draw from any 
finite-dimensional distribution of the target SDE by rejection sampling with 
proposals from the Wiener measure. In [8] it was noticed that the transi- 
tion density of the target diffusion can be written in terms of the acceptance 
probability of EA, thus, a simple pointwise estimator of L{9) for any € is 
readily available. The algorithm is termed the Acceptance Method (AM) in 
[8]. In that paper a simultaneous version of the algorithm is also presented, 
but for the case of a specific family of diffusions where its development (from 
pointwise to function estimator) is rather straightforward. In the current pa- 
per SAM is applied to a considerably larger family of diffusions and its con- 
struction will involve novel couplings of Brownian motion paths. Also, the 
consistency properties of the method (conditional and unconditional, i.e., 
for fixed and of increasing number data points respectively) are examined 
here for the first time. 

The applicability of SAM (and AM) is attached to that of EA, whose 
latest development [7, 8] covers the class of diffusion processes determined by 
the set of conditions (Co)-(C3) on the drift and diffusion coefficient given in 
Section 2. Though not universally applicable, SAM can be applied to several 
diffusion models for which the likelihood function is intractable, for example. 
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to the class of diffusion processes generated by one-to-one transformations 
of the hnear diffusion with multiphcative noise (see, e.g., page 119 of [29]): 

(5) dVs = {6i + e2Vs) ds + (^3 + O^Vs) dBs, 

with state space V = (max{— ^3/^4, —6i/92\-,oo). We will discuss the poten- 
tial of SAM for more general classes of diffusions. Relevant to this direction 
are undergoing developments in EA [6]. 

The paper is organized as follows. In Section 2 we set up the basic nota- 
tion and provide the random function •); we state its properties in The- 
orem 1. In Section 3 we establish a.s. uniform convergence of the estimators 
£^(-) to £n(')i ^ N^oo, and discuss some important consequences. In 
Section 4 we illustrate our method by applying SAM to an example SDE. In 
Section 5 we allow n — > 00 and present the theory suggesting computational 
optimality for the choice N = 0(n^/^). In Section 6 we prove Theorem 1. 
In Section 7 we summarize and discuss further ideas. The paper is supple- 
mented with a brief Appendix. 

2. Basic notation and statement of the main result. In this section we 
set up the basic notation, construct explicitly the estimator L(H, •) of L(-) = 
Pt{v, w] •), and state the main, from the applied perspective, result of the pa- 
per, which is proven later in Section 6. Throughout the rest of the paper we 
assume that is a compact subset of R'^ which contains the unknown MLE 
On - Generally, random variables will be written in capital letters; typewriter 
style will be used to emphasize that the distribution of the variable is inde- 
pendent of 0. 

We define 




ueV, 



to be any fixed anti-derivative of l/(7{-;6). For arbitrary u, let 

a{r] ^{u,0);d) Jq 

where r]~^ is the inverse of rj. Assuming that r]{-,9) is twice continuously dif- 
ferentiable, Ito's lemma shows that Vg > rj{Vs;0) =: Xg is the unique trans- 
formation (up to a change of sign) that maps y to a process of unit diffusion 
coefficient; a{-; 6) is the drift of the transformed process X. Throughout the 
paper we assume that the following conditions hold for all € 0: 

(Co) '?(•; ^) is twice continuously differentiable; the law of Xg = rj{Vs,9) has 
a density w.r.t. the Wiener measure provided by Girsanov's theorem 

(13); 

(Ci) a{-;9) is continuously differentiable; 
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(C2) (a^ + a'){-;9) is bounded below; 

(C3) (a^ + a'){-;9) is bounded above on (z, 00) for all 2; G R. 

Remark 1. The essence of Theorem 1 below is to devise an unbiased 
estimator of the transition density of the unit diffusion coefficient process 
Xg =r]{Vs,9) when its drift a{-;6) satisfies conditions (Ci)-(C3) above. 
Thus, the case when (a^ + «')(•; ^) happens to be bounded above on {—00, z) 
is covered by symmetry. One then needs only to consider instead the (unit 
diffusion coefficient) process —77(1/5,^); its drift will then satisfy (Ci)-(C3). 

We define 

m = mf + <^(., 9) = (J^±^^ _ 1^0) > 0, 

and r{u, 9) = su^^(z(^u,oo) 4>{z, ^) < ^■ 

The random element in the estimator L(H, ^) is H = (E, Z, N), where E ~ 
Ex(l), Z~A/'(0, 1), * is a homogeneous Poisson process on [0,t] with time- 
ordered points Yj, 1 < j < A, of number A Poisson(Ai), where the intensity 
A is specified below, and conditionally on A, the 3 x A-matrix N = {N.^. , 1 < 
* < 3, 1 < j < A} consists of independent standard Gaussian variables. On 
the event {A = 0}, * and N contain no elements. The intensity A depends on 
E and the observed data, v,w, in the following way: we define the functions 

(7) x = x{6):=r]{v,6), y = y{9) ■= r]{w,6), 

(8) m = m(E, 9) = {y + x - sj2tE + {y - x)2)/2, 

and take A to be any real not less than supggQ r(m, 9). The following theorem 
gives the random function L(H, •) as a composition of easily computable 
functions. We define Aft{u) := e-"'/^^*) u G R, t > 0. 

Theorem 1. Under assumptions (Co)-(C3) and the following regularity 
conditions: 

(Cntl) (i) a{-;-), a' {■;■), and (ii) A{-,-) are continuous on R x Q , 
(Cnt2) (i) r]{u,-), (ii) r]'{u,-) are continuous for all u £ V , (Hi) l{-) is con- 
tinuous, 

the random function L{E,-) defined below is such that: 

(i) for any9€G, E[L{E,9)] = L{9), 

(ii) w.p.l 9 L{'E,9) is continuous, 

(iii) there exists a constant M such that w.p.l \L{E, 9)\ < M for all 9 £ @. 
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L(H, e) = \7]'{w, 6)\J\ftiy - x) exp{A(y, 9) - A{x, 9) - l{9)t} x a(H, 9); 

a{Z,9) = j2m^\{[l-^{x.,,9)/\]- 
i=i j=i 



X„ =m + 



\ 



t 3 \ I ^ \ ( ^ 

+ E 7., J + E 7., J + E ^3. 7,,. 
^ i=\ ) \i=\ I \i=\ / 



/?,^ = (x - m)^^^I[Yj < n\ + (y - m)^5^I[Yj > r^]; 

T% t-Ti 

, _ (r.-Y,)^(Y,-Y,_0 ^ 
(r.-Y0(r.-Y,_O '^^^ " 

(t-Yj)^(Y;-Y;_i VtQ ^^^^ 

Ti = t[l-\ g] , T2 = t[l^ g 

V X — m / V X — m 

_ tgE + 2{x-mY 

(l + <7)[tE + 2(x-m)2]' ^^2 = 1 -Pi; 



5 = l + Z^/E - y'2Z2/E + Z4/E2 >0. 

Remark 2. Notice that Xij: Pij^ liju "^j; Pi functions of 9, the 

random element H and the data. Also, H depends on the data points v,w 
through the Poisson rate A. 

Corollary 1. Under the assumptions of Theorem 1, the likelihood term 
L(-) is continuous. 

Proof. Consider some ^ G and a sequence {9j} such that 9j — > 9. 
Then, lim^^oo L{9j) = limj^ooE[i(H, 9j)] = E[limj^oo L{E, Oj)] = L{9), where 
the first equality holds due to the unbiasedness stated in result (i) of Theo- 
rem 1, and the others due to results (ii), (iii) and the bounded convergence 
theorem. □ 



3. Uniform convergence of likelihood estimator. The Monte Carlo av- 
erage Lf [9) in (3) is an unbiased estimator of the likelihood factor Li{9) 
for any 9 E @, thus, Kolmogorov's Strong Law of Large Numbers (SLLN) 
implies that w.p.l, £^ (9) £n(^) as ^ oo. It is known, however, that 
this pointwise convergence is not strong enough to guarantee convergence 
of the maximizers 9^ of 2,^ {■) to the MLE Similarly, a stronger form of 
convergence is needed to ensure that several other interesting features of the 
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Monte Carlo functional averages, such as level sets, integrals over subsets 
and profile likelihoods, will converge to the corresponding features of the 
likelihood function. A sufficient condition is that the convergence is uniform 
in 6: 

(9) w.p.l, lim sup|i:^(0) =0. 

Questions which lead to essentially equivalent problems to uniform con- 
vergence of functional averages occur in theoretical statistics (see, e.g., [11, 
38, 39] for the Glivenko-Cantelli problem and [13, 40] for the consistency of 
the maximum likelihood estimator), dynamical systems and ergodic theory, 
stochastic optimization, econometrics [2] and Monte Carlo methods [22]; see 
Chapter 1 of [33] for a review. The general probabilistic framework in which 
convergence of functional averages is most naturally addressed is probability 
on Banach spaces. 

We prove (9) using the following fundamental theorem about SLLN for 
random elements in an arbitrary separable Banach space. This result first 
appeared in [30]; see also [5] and [24]. We recall (see, e.g., [24]) that for an 
arbitrary random element X m. a Banach space (X, || • ||), \\X\\ denotes its 
norm, and when E[||^||] < oo, the expectation E[A'] is defined as the unique 
element ^i€X such that T{fj,) = K[T{X)] for every linear function T : X ^ R 
in the dual space of X. 



Theorem 2. Let {X, 
element in X, such that 



be a separable Banach space and X a random 



^[\\X\\]<oo- ¥.[X] = 0. 
If Xi,X2, ■■ ■ are independent copies of X , then 

N 



W.p.l, 



lim 



-yx. 



The uniform convergence in (9) follows easily from the following corollary. 



Corollary 2. Let L[=.,d) be the random function constructed in The- 
orem 1, and "E} ^ . . . , independent copies ofE. Then 



w.p.l, 



lim sup 



N 



yL{-^,e)-L{e) 
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Proof. Take (X, || • ||) to be the space of continuous real functions on the 
compact set Q equipped with the sup-norm such that, for any element / G 3£, 
11/11 = supgge 1/(6*)! . This is well known to be a separable Banach space. 
Theorem 1 has shown that L(E,-) takes values in X, and E[||L(E!, •)||] < oo. 
Corollary 1 has shown that L(-) G X, and by uniqueness of expectation, 
E[L(5, •) - L(-)] = 0. Applying Theorem 2 to L(H, •) - L(-) yields the result. 
□ 

The compactness of and (9) imply the following result, which validates 
our Monte Carlo maximum likelihood approach. 

Corollary 3. If On is the unique element o/ argmax^gg and 
{0n}N any sequence of maximizers o/{£^ (Oliv; thenlim.N^c<yOn = w.p.l. 

Note that the same result holds for a sequence of so-called e^^'-maximizers, 
that is, a sequence O^e such that — 6^^\ < e'^^ , where e^^' — > as ^ oo. 
This extension allows for numerically efficient implementations of SAM, as 
in Section 4. 

Certainly, consistency can also be established under the classical approach 
for showing convergence of the MLE from i.i.d. data (because of the indepen- 
dence of the Monte Carlo samples in the case of SAM) to the true parameter 
value; see, for example, [19]. It is our impression though that a generalized 
SLLN provides a natural approach for comprehending strong convergence 
in enlarged spaces. 

Note that a.s. continuity of the random function is not a necessary con- 
dition for uniform convergence. For instance, concavity of the functional 
averages and their expectation would suffice; see [23]. Uniform convergence 
might hold even when H^,H^, . . . , form a stationary (not necessarily ergodic) 
stochastic process; see, for example, [13, 15, 22, 33]. On the other hand, 
consistency results as in Corollary 3 could be obtained by establishing epi- 
graphical/hypographical convergence (see, e.g., [4] for relevant results from 
convex analysis) which is weaker than uniform, although closely related, 
and requires only semi-continuity properties of L(H, •); see [22] and [13] for 
investigations in the context of functional averages. In our case, however, 
concavity does not necessarily hold, but continuity is rather easily estab- 
lished. 

Having shown uniform convergence, we can resort to Theorems 5 and 6 of 
[22] to prove that the Monte Carlo profile log-likelihood and the Monte Carlo 
level sets converge to the corresponding features of the likelihood function 
as the number of Monte Carlo samples increases. 
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3.1. Asymptotic normality. Having established strong consistency of 9^ 
as an estimator of the unknown MLE 0„ as N ^ oo, we can proceed, under 
some additional conditions, to prove asymptotic normality for the rescaled 
sequence N^/^{e{^ - On). 

We will establish this result by appealing to a general theorem stated 
below, which concerns the asymptotic normality of the maximizer of an un- 
biased simultaneous estimator of the likelihood, as the Monte Carlo sample 
size — > oo . We first state the result in a general context and then discuss 
when the conditions it requires are satisfied in our context. We will need the 
log-likelihood function 

n 

ln{0):=Y.\ogU{e) 

i=l 

and its estimate 

n N 

ie):=Y: log L fie); Lf{e) = y: UiEi,e)/N. 
i=l j=l 

The random elements h| , . . . , are independent copies of, say, Hj , which is 

used for the estimation of the ith likelihood term Li{6) = PAt^ (Vt._i , Vf. ; ^). 

The random elements are independent over the data point index i = 1, . . . , n. 
v c 

We use — >, ^ to denote convergence in probability and distribution respec- 
tively. All derivatives in the sequel are w.r.t. the parameter argument. 

Theorem 3. Assume the following: 

(a) On is unique, in the (assumed nonempty) interior of Q. 

(b) e^^On asN^oo. 

(c) L{9) =E[L(H,0)] can he differentiated twice under the expectation 
sign. 

(d) There is convergence in distribution 

for some covariance matrix An. 

(e) Bn = —V^inidn) is positivc definite. 

(f) V^i^ (6) is bounded in probability uniformly in a neighborhood of On- 
It is then true that as N ^ oo, 



and 

N^i\e^ - en)^M{o,B-'AnB~'). 
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Proof. This is a known result; see, for instance, [19, 22]. Briefly, for 
the first result note that (b) and the uniform bound in (f) imply, from 
the mean value theorem, that V^£^{e^) - V'^i^iOn) ^ in probability. 
Prom (c) and SLLN we get that a.s. V'^i^ {On) V^in{6n)- Por the second 
result one takes a first order Taylor expansion with integral remainder for 
VC(^n ) - ^C0n) and rescales by iVVs. □ 

The most restrictive condition of this theorem in the context of SAM 
is the differentiability of 6 ^ L(E,6). The formula in Theorem 1 suggests 
that L{E,-) is typically nondifferentiable at {9 £ @:Ti{6) = Yj,i = 1,2, j = 
1, . . . , A}. A simple case where this is guaranteed to be an empty set is when 
the diffusion coefficient a does not depend on 9. If cr{u; 9) = ct{u), then also 
r/(ii, 9) = r]{u), and 9 is involved in 9) only through the second argument 
of (t>{u, 9). Differentiability for L(H, •) can now be implied by straightforward 
conditions on 9 ^ a{u;9). Thus, Theorem 3 applies directly to SAM under 
the assumption that (t{u; 9) = a{u). Note, however, that we have experimen- 
tally verified consistency of the same rate 0{N^^'^) even for the general case; 
see, for example, the numerical example in the next section. The points of 
nondifferentiability are of a.s. finite number, and their presence does not 
seem to effect the result (of local character anyway) of Theorem 3. 

We exploit the independence of the transition density estimators Lj(Hj,0) 
over i = 1, . . . ,n to establish condition (d) and identify An- Notice that 



ivV2v£^(4) = v(Ari/2f 



i=l ^ ^ Lf{(^n) Li{On 



Proposition 1. Let {Xj,Yj} he an i.i.d. sequence of vectors Xj and 
positive scalars Yj with expected values Hx and fiy respectively, and running 
averages Xn = EjLi /N, Yn = Y.f=i Yj /N. When N^oo, then 

N^^f^ - !^)^^M{0,YaT{f,yX, - Y^^,x)). 

\Yn fly J fif, 

Proof. It follows directly from Slutsky's theorems; see [19]. □ 

Using this proposition, we find that 

^ _ ^ Var(L,(g„)VL ,(H,, 9^) - L,(H,, 4)VL^(4)) 



i=l 



Li{9r, 



^Varfv^^i^=^ 

i=l \ Li{On) 
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An can be estimated via the simulated sequences {Lj(:;:- , 0^), VLj(:i^, 
Note that A n increases only linearly with n. 

4. Numerical illustration. We applied SAM to the logistic growth SDE 
[7]: 

(10) dVs = 5Vs{l-c-^V,)ds + aVsdBs, eV=(0,oo), 

which is used to model the evolution of a population in an environment of 
capacity c; 5 is the rate of growth per individual and cr > a noise param- 
eter. It is known (see, e.g., page 123 of [29]) that (10) is the inverse of the 
linear SDE with multiplicative noise (5) after making the correspondence 
{01,62,03,64) = {5/c, cj^ — 6, 0, —cr). Here, we take 6 = {6, c, a), and 

r]{u, 6) = - log(u) /a, a{u; 0)=a/2- 6 /a + V(c7c)e-'"", 

1(9) = a^S - 6/2, r{u, 0) = [{a^ + a'){u; 0)/2 - l{0)] V [6^/{2a% 

Note that, following Remark 1, r]{u, 6) is defined here as the negative of the 
transformation (6). If G = x [QiCu] x [o";,(Tm], then for any given pair 

of successive data points v,w of time increment t, 

A = 7^2 X [(e'/'/Q - 1)' V 1]; q := log{vw) + ^2talE + log'{w/v). 

It is easy to verify that all conditions of Theorem 1 hold. 

In Figure 1 we demonstrate numerically the consistency of the Monte 
Carlo MLE as an estimator of the unknown MLE. We applied SAM to a 
dataset of size n = 1000 under the specifications Vq = 700, 6*0 = (0.1, 1000, 0.1) 
and Ati = 1. (The dataset was simulated using the Exact Algorithm of [7].) 
We chose 6 = [0.03,0.18] x [850,1200] x [0.09,0.12], outside which a pre- 
liminary investigation showed that the likelihood is of negligible value. The 
maximizers 6^ of 2,^ (') foi' the various values of in Table 1 were found 
numerically with the downhill simplex method (see, e.g.. Section 10.4 of [34]) 
up to some precision error e*^' which decreased with increasing N. The ini- 
tial search point for each maximization was the output of the previous one; 
for A = 1 the initial point was (0.05,1150,0.115). Computer implementa- 
tion in C on a Pentium IV 2.6 GHz processor yielded (in 14 minutes) the 
sequence of e'^^-maximizers shown in Table 1. The simulated A was on aver- 
age, over all consecutive data points and Monte Carlo iterations, 2.004 and 
never exceeded 11. 

In Table 2 we investigate the mean of the asymptotic distribution of the 
scaled variable N^/'^{On — On) theoretically demonstrated in Theorem 3. Note 
that in this context the diffusion coefficient depends on unknown parame- 
ters, so ^ 1-^ L(E,0) will exhibit points of nondifferentiability with positive 
probability. Yet, the numerical study suggests agreement with the conclu- 
sions of the theorem. 
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Table 1 

Consistency, A'' — > oo; The e'-"'' -maximizers of the estimated likelihood 2,^ (■) for a dataset 

of size n = 1000 from the logistic growth model simulated under the parameter values 
6o — (0.1, 1000,0.1). In parenthesis the standard errors found by inverting —V^in{dn,e)- 



N 


"rt,e 






1 


0.1063 (0.01493) 


1010.0 (30.25) 


0.10051 (0.002354) 


2 


0.1105 (0.01539) 


1009.8 (29.67) 


0.10053 (0.002352) 


5 


0.1115 (0.01562) 


1010.4 (29.31) 


0.10057 (0.002369) 


10 


0.1103 (0.01559) 


1012.3 (29.84) 


0.10060 (0.002364) 


20 


0.1097 (0.01558) 


1012.3 (29.99) 


0.10057 (0.002372) 


50 


0.1100 (0.01563) 


1012.9 (30.02) 


0.10058 (0.002374) 


100 


0.1099 (0.01563) 


1014.2 (30.07) 


0.10058 (0.002371) 


200 


0.1095 (0.01559) 


1014.4 (30.19) 


0.10057 (0.002368) 


300 


0.1096 (0.01560) 


1014.4 (30.19) 


0.10057 (0.002367) 


Large N{= 10*) 


0.1096 (0.01561) 


1014.5 (30.18) 


0.10057 (0.002366) 



5. Unconditional asymptotics, n — > oo. We find the optimal choice of 
Monte Carlo iterations N = Nn asymptotically as n — > oo. In accordance, 
we will now write for the sequence of Monte Carlo maximizers (4) and 
Lf" (9) for the unbiased estimators of Li{9). The data points are now treated 
as random variables. The analysis follows a rather intuitive approach. The 
objective is to provide a rule for the selection of N but, at the same time, 
avoid the mathematical rigor that would take up the space of a full paper; 
see, for example, [27] or [12]. So, we will state the main result in Theorem 4 
under reasonable ergodicity assumptions on the diffusion dynamics without 
insisting on technical details. Note that the precise formula for the estimator 
L{'E,0) of L{6) is not essential for this section: given that Assumptions 1- 
4 below are satisfied. Theorem 4 holds for any given unbiased estimator 
of L{6). As in Section 3.1, the asymptotic normality result in Theorem 4 
requires differentiability of i— > L{E, 6) and will apply directly to SAM under 
the known diffusion coefficient assumption, a{u;9) =a{u). 



Table 2 

Asymptotic Unbiasedness, N ^ oo: The sample means from 1000 realizations of 
N^^^{9^ — On) for various N . In parenthesis the corresponding standard errors. The data 
were of size n = 250; Oo = (0.1, 1000,0.1). The MLE ^„ was found using large N (= lO"). 







N''\5^ - 






.) 


N = 


25 


-129e-5 (30e-5) 


0.825 (0.320) 


-2.28e-5 (0.79e- 


-5) 


N = 


50 


-81e-5 (30e-5) 


0.414 (0.321) 


-0.47e-5 (0.78e- 


-5) 


N = 


100 


-46e-5 (30e-5) 


0.092 (0.317) 


1.40e-5 (0.81e- 


-5) 
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We assume that the diffusion is ergodic, with vr denoting its invariant 
density corresponding to the correct parameter value and that the data 
are equidistant, that is, ti = iA for some A > 0. We define the family of 
mappings which satisfy a Law of Large Numbers criterion: 

9 = |/I/:R'^R.,E fiyii-i)A,ViA)/n ^ E,[/(yo, ^a)] | , 

where K^^l] is expectation in stationarity, that is, {Vq, Va) ~ TT{dx)pAix, dy; 9q). 

We follow the approach of [27]. In that paper the discretization increment 
of a diffusion approximation is adjusted to the datasize n — > oo. In our case, 
had it been possible to construct an unbiased estimator of the Zo^-transition 
density, then the Monte Carlo error would be averaged out for increasing 
n and -y/n-consistency of 9^" (as an estimator of ^o) would follow even for 
fixed N; see Theorem 2 of [27]. We will now allow N = Nn — > oo as n — > oo 
and identify the magnitude of the log-bias through an error expansion. For 
simplicity, we set 

Ci = Ci,N := VlogLf (0o) - VlogLi(0o), 

and J^i = o"(V(j_i)A) ^ia), for z = 1, . . . , n. Recall that d is the dimensionality 
of the parameter vector. 

Assumption 1. There exist tl^ -.K"^ ^ BJ^ , g : R'^^'^, /iiR^^^R 

with scalar components in 9 such that 

E[Ci,N I - m{i~l)A,ViA)^ + nCi,NCl,N I -^d " 5(^(i-l)A , V^^a) ^ 



N 

<hiV^i_,)A,ViA)o{l/N). 



The assumption is not on the expansion itself, but on the regularity of the 
coefficients ip,g,h. Indeed, for given V(j_i)A,ViA5 the 0(l/A'")-bias follows 
from a second-order Taylor expansion for Ci,N- Analytically, for sequence 
{Xj,Yj} as in Proposition 1, one can write 

^ - — = —{Xn - fix) - ^{Yn - fly) - \{Xn - flx){YN - fJ-y) 
IN fJ-y fJ-y f-y fJ-y 



+ —{YN-fly) + R{XN,YN,flx,fJ.y), 

for some lower order residual R{XN,Yiy, fi^, fiy)- Making the correspondence 
Xn ^ VL^{9o), Yn ^ Lf{9o), assuming that E[Li{Ei,9)] = Li{9) can be 
differentiated under the integral sign, and taking expectations conditionally 
on the data, ?/'(^(i-i)Ai ^a) in Assumption 1 is analytically identified as 

^ -Cov(VL,(H„0o),^i(Hi,eo)) + ^fe7rvVar(Li(Hi,0o)). 
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The corresponding expressions we can obtain for g, h are much more com- 
phcated. 



Assumption 2. The matrix A{9) := E7r[-V^ logpA(Vb, Va; 6*)] is posi- 
tive definite for any 9 £ Q and 



n 

uniformly for 9 in a neighborhood of 9q. 



Under the standard maximum hkelihood assumption for convergence in 
probabihty of 'Sl'^ln{d)/n to A{9) uniformly in a neighborhood of it re- 
mains to explain the same mode of convergence for {V^^^" {9) — V^^n(0)}/n 
toward 0. We define 

6 = C.,7V = V^logLf (0) - v2logL,(^). 

From Lemma 9 of [21], used also in [27], the following conditions imply the 
required convergence for fixed 9: 

n ' 

Similarly to Assumption 1, a Taylor expansion can provide estimates: 

\E[C,,N I ^i]| + nCiM^N I H\ < Re{V{^-l)A,V^A)0{l/N). 

Assuming that Rg and since Nn oo, these estimates imply (11). It is 
technically much harder to illustrate convergence uniformly in a neighbor- 
hood of ^0- Following Theorem 2 in [27], it suffices (assuming stationarity) 
to obtain a bound on the L^'^"'"^-norm of Ci,Ar(^) and a Lipschitz condition 
for 9 I— > Ci,Ar(0) uniformly in N, again in the L^'^'^^-norm. We avoid further 
details. 

Assumption 3. The following two sequences converge in probability to 

0: 

n ' n 



Note that the summands have zero expectation. This can therefore be 
interpreted as an ergodicity assumption on the diffusion dynamics upon the 
consideration of the enlarged cj-algebra from the Monte Carlo scheme. 
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Assumption 4. The following weak convergence holds: 

where 

This is a standard ergodicity assumption for maximum likelihood infer- 
ence for discretely observed diffusions. Of course, V = A{0q), assuming ex- 
changeability of differentiation and integration at ^^^[V^paI^j ^A! ^o)]- 

Theorem 4. Let N = Nn grow to infinity with oo. Then 

Also, if 9q is in the interior of Q and Assumptions 1-4 hold, then: 

(i) if lim„^oo V^/Nn = then ^{6^- - 9o) ^ AA(0, A~^), 

(ii) i/lim„_,oo\/^/iVn = cG (0,cx3) then ^/E{d^'- - Oq) ^J\f{cn,A~^), 

(iii) if lim„^oo y/n/Nn = oo then NniOn" - Oq) ^ 
forA = A{9o) and 12 = A-^E^['l|){VQ,VA)]■ 
PROOF. Consistency of the Monte Carlo MLE can be proved in line 

with [12]. We proceed directly at the proof of the asymptotic results for the 
various scalings of Nn. We use a Taylor expansion as in the classical MLE 
theory. The additional complexity is due to the presence of the Monte Carlo 
scheme. The basic equation is the following: 

(12) (9^) = vC{9o) + C VH^ {00 + p{9^ - 9o)) dp x {9^ - 9o), 

Jo 

where the term Vi^{9^) = can be ignored. From the consistency of 9^" 
and Assumption 2, 

- V^e- {9o + pi9 ^--9o))dp V ^ 



n 

Consider now the remaining gradient term. We rewrite 

Vi^ i9o) = {V^ i9o) - V£ni9o)} + V4(^o). 

Assumption 4 controls V£n(6'o). For the last term {Vi^ {9o) - V£„(6'o)} we 
employ a martingale decomposition as in [27]. We rewrite 

n n 
1=1 i=l 
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Assumption 1 controls the extreme right term. Under Assumptions 1 and 3, 

n ' 
so the martingale CLT in Theorem 3.2 of [26] implies 

Er=i{Ci,Af„ - HCi,N„ I ^i]} c Q 

For (i) and (ii), multiply both sides of (12) with whence, from As- 

sumption 1, the log-bias term X^i^i IE[Cj,JV„ I ^i] will converge in prob- 

ability either to if n^^'^/Nn — > or to cfi if n^^'^/Nn — > c. For (iii), multiply 
both sides of (12) with Nn/n. □ 

The theorem shows that, as long as = o(n^/^), the Monte Carlo MLE 
has the same asymptotic behavior with the unknown MLE. 

6. Proof of Theorem 1. The notation in this section follows the defini- 
tions in Section 2. We prove the result in several stages. 

6.1. Diffusion transformation and densities. We consider the modified 
process Xg = rj{Vs,6)- Condition (Co) allows the application of Ito's lemma 
to show that Xg is the solution of the unit diffusion coefficient SDE: 

dXg = a{Xs;e)ds + dBs. 

Let pt{-, •; 6*) be the transition density of X defined analogously to (2). Recall 
that we have defined x = x{6) = r]{v,9) and y = y{9) = r]{w,0). A standard 
change-of-variables argument yields 

HG) = Pt{v,w-e) = \r]'{w,9)\ ■pt{x,y-e). 

Let Qe be the law of the paths of X on [0, t] conditioned to begin at Xq = x 
and finish at Xt = y. We call such a path "a bridge from {0,x) to {t,y)." 
Let We be the distribution of the Brownian bridges (BBs) from (0,a;) to 
(t,y); a random process distributed according to We will be denoted by W. 
Both Qe and We apply on the space C of continuous mappings from [0, t] 
to R equipped with the corresponding cylinder cr-algebra; we denote by u! 
a typical element of C. We have assumed in (Cq) absolute continuity of the 
law of (unconditional) paths of X w.r.t. the Wiener measure with density 
given by Girsanov's theorem (see, e.g., [31]): 

(13) exp|y a{uJs',d) diJs — J a^iuJs'jO) ds^ . 
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Bayes' theorem and an application of Ito's formula gives the foUowing ex- 
pression for the corresponding bridge density: 

(OJ) — -1'-. ^^ 



■ exp|A(y,e) - A{x,e) - ^ ^(a^ + a')(a;.;0) dsj. 



(iW/ ' Ptix,y;e) 

Taking expectations at both sides w.r.t. Wg, using (C2) and re-arranging, 
we get that 

pt{x, y- 6) =Mt{y- x) exp{A(y, 6) - A{x, 6) - l{9)t} x a{9), 

where 



a{e) = E ■ 



exp| — y (j){ujs,9) ds 



< 1. 



A heuristic description of the remainder of the proof is as follows. We 
initially derive an unbiased estimator a{3,6) of a{9), where the distribu- 
tion of the random element H depends on 6. We then construct jointly the 
family G G} by expressing H = /(H,X,0) for H given in Theorem 1, X 

some other random element also independent of 9 and / an appropriately 
specified function. Finally, we show that a{'E,9) defined in Theorem 1 is 
given as a{E,9) =K[a{f{E,X,9),9) \ H], where X is integrated out to ensure 
a.s. continuity of the mapping 9 a{E,9). 

6.2. Connection with exact diffusion bridge simulation. In [8] it is no- 
ticed that a{9) coincides with the acceptance probability of a rejection sam- 
pling algorithm for the simulation of paths with law Qg. The algorithm, de- 
veloped in [7] and termed the Exact Algorithm (EA), proposes paths from 
Wg and accepts them according to the density ratio 

(14) ^(^)ocexp|-yV(^s,^)rfs|<l. 

Let m = inf{ws,s G [0,i]} be the minimum of to. Condition (C3) implies 
that (j){ujs,9) /r[m,9) < 1, for all s E [0,t]. When u; is a realization of the 
Brownian bridge W ~ Wg , the distribution of its minimum is given in terms 
of the Rayleigh distribution (see, e.g., [36]) and can be simulated precisely 
according to (8) using the exponential random variable E. Notice that if $ 
is a Poisson process of intensity r{m,9) on [0,t] x [0, 1] and N the number 
of the points of $ lying below the curve s 1— > (j){uJs, 9) /r{m, 9), then 

P[iV = I w] = exp|-y 0(w,,6')ds|. 

Thus, in [7] the following simulation algorithm is suggested. Let ^ be the 
projection of $ on the time-axis with time-ordered points Yj,l <j < A, 
where A Poisson(r(?7i, ^)t). 
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The Exact Algorithm (EA). 

1. Simulate E ~ Ex(l), set m = (y + x - ^2tE+ {y - xf) /2. 

2. Simulate 

3. Simulate {WVj , 1 < i < A}, so that its minimum is m. 

4. With probability 1 - 11^=1 [1 - 0(Wy.,6')/r(m,0)] goto Step 1. 

5. Output all information about W . 

Step 3 requires the simulation of the skeleton {Wy^. ,1 < j < A} with a 
given minimum. This can be achieved (see Section 6.4 below) using existing 
theory on the decomposition of the Brownian path at its minimum; see, for 
example, Proposition 2 of [3]. Therefore, a pointwise unbiased estimator of 
a{9) is readily available: 

A 

(15) a(H, e) = l[[l- <P{Wy, , e)/r{m, 0)] , 

j=i 

where H = (E, {VFy. , 1 < j < A}). This is the Acceptance Method proposed 
in [8]. 

6.3. Coupling of the Poisson processes. We define the joint structure of 
the collection of random variables {^;9 € Q} using the thinning property of 
the Poisson process (see, e.g.. Section 5 of [28]). Let * be a Poisson process of 
rate A > supg^Qr{m{E,9),6) on the interval [0,t] with time-ordered points 
Yj, 1 < j < A and U = (Ui, . . . ,Ua) a vector of i.i.d. variables, Ui ~ Un[0, 1]. 
We set ^ = {Yj e * :Uj < r(m, 9)/X, 1 < j < A}. Then the right-hand side of 

(15) rewrites as 

A 

I[Uj- < rim, e)/X] ■ cl){W^^ , e)/r{m, 9)}. 

i=i 

After integrating out U, we have the following pointwise unbiased estimator 
oia{9): 

A 

(16) a{E,9) = l[[l-cl>iW,^,9)/X], 

where now E = (E, *, {W^. ,1 <j < A}). 

6.4. Coupling of the proposed paths. For a path to with given minimum 
m, let T = {s £ [0, t]:ujs = m} be the instance when m is attained. As shown 
in the Appendix, when u; is a realization of W, the time instance of the 
minimum can be simulated conditionally on m as follows: 



(17) 



T = I[Y < pi]Ti > pi]t2, 



MONTE CARLO MAXIMUM LIKELIHOOD FOR DIFFUSIONS 



19 



for V Un[0,l], where we recall that Pi,ti,T2 are given in terms of E, Z 
and 6. A Brownian bridge W with a given minimum m attained at time 
r can be constructed in terms of two Bessel bridges which can in turn be 
obtained through six independent standard [from (0,0) to (1,0)] EEs; see 
[3]. In detail, the theory allows us to construct W at the time instances 
Yj, 1 < J < A, in the following way: 




m)(r 



^3/2 



+ Wi,Y^/, +W: 



+ W: 



(18) Wy, =m+< 



{y-m){Yj - 
{t - r)3/2 



4,(Y, -r)/(i-r) 



+ W5,(Y,-r)/(i-r)+W6,(Y, 



-r)/(t-r)2 



3,Y,/r' 

if Y,- < T, 



1/2 



if Y,' > r. 



where W = {(Wj^^, < s < 1), 1 < i < 6} is the required collection of standard 
EEs. We define VFjjj as W^^ in (18) under the specification r = r^, for i = 1, 2 
and 1 < J < A. Under this convention, substituting (17), (18) into (16) and 
integrating out V, we rewrite the function we wish to estimate as 



a{e)=Em 



n[l-</'(W^Y,,^)/A]|E,Z,W,* 
■i=i 



= E ^1 n [1 - <^(^^Y, , e)/\] xp\= ai{9) + a2{e), 
1=1 [j=i ) 

where ai{6), a2{0) are defined in the obvious way. In the sequel, we construct 
estimators aj(H, 6) of ai{6),i = 1, 2, and take a(H, 9) = ai(E, 6) + a2(H, 9). 

For each Yj, we need to simulate the bridges Wi, W2, W3 at the instance Yj/rj 
or simulate W4,W5,W6 at the instance (Yj — rj)/(t — Tj) if Yj < ti or Yj > Tj 
respectively, 1 < j < A; see also Figure 1 for a graphical illustration. Thus, in 
total we need 3 x A simulations, which can be done using the same number 
of Gaussian random variables as we now show. We can simulate a standard 
EE (W^; < s < 1) at some time instances < si < • • • < < 1 using the 
formula (derived using standard properties of EE) 



(19) 



1=1 



(l-Si_i)(l-sO 



1 < i < 



where Ni,N2, ■ ■ ■ ,N(i are i.i.d. with A'^i ~AA(0,1). Consider now the col- 
lection N = {Nyfcj, 1 < A: < 3, 1 < j < A} of independent AA(0, 1) variables. We 
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Wfc 

n ^ 










r ' 




^ ■ ■ % 
















i-^ \ 1 


T^- +' if 1 


!^ ■ 








t-Tj '-tf 



Fig. 1. Scheme for simulating and \lk+3,k — 1,2,3, at the required discrete time in- 
stances; we have defined hi = X/j-i ^I'^i — '''»]■ 



will use {Nfcj, 1 < J < A} for the generation of the required instances of the 
bridges and 1 < A; < 3; see Figure 1. Following (19), we obtain 



k + 3,{Yj-Ti)/{t-Ti) 



for 1 < i < A, A; = 1, 2, 3, and for both i = 1, 2. Notice that there are alterna- 
tive choices for the realization of the required locations of the standard BBs. 
However, the above formula has been carefully developed to ensure continu- 
ity in the parameter 9. These expressions allow us now to write the proposed 
paths as deterministic functions of E, Z and N. In particular, substituting 
the above expressions into (18) and setting Xij equal to the right-hand side 
of (18) for r = n, with i = 1, 2, 1 < j < A, we obtain the final estimator: 

A 

ai(E,e)=pixl[[l-<p{x,^,e)/\], i = l,2. 

J=l 

6.5. A.s. continuity of the random function. (Cnt2)(i) implies that 74^7 
are continuous functions of 6 for all i,j,l (even at values of 9 such that Yj = ti 
or Y; = Ti); so, is continuous in 9. (Cntl)(i), (Cnt2)(iii) suggest that (/>(-, •) 
is continuous on R x 0, thus, (l){Xij-,S) and a{E,9) are continuous in 9. The 
continuity of L(E,9) follows from assumptions (Cntl)(ii) and (Cnt2). 



1=1 



/ (r,-Y,)2(Y,-Yz:i) 

Ti{Ti-Yl){Ti-Yl_i)' 



1=1 



> (r^-Yj)2(Y;-Y;_i VrQ 

{t-Ti){Ti-Xl){Ti-Yl_iyTi) 

> Ti, 



II[Yi >Ti], 
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6.6. A.s. boundedness of the random function. Note that a{=.,0) < 1 a.s., 
thus, w.p.l \L{E,e)\ < M for all 6* G 9, where 

M = sup{\7]'{w, e)\Nt{y - x) exp{A(y, 6) - A{x, 9) - l{e)t}} < oo. 
6»ee 

M is finite as the supremum of a continuous function over a compact set. 

7. Discussion. In this paper we have introduced a new method for like- 
lihood inference for discretely observed diffusions. The method is compu- 
tationally efficient and simple to implement, and expands significantly the 
family of diffusion models for which routine maximum likelihood calcula- 
tions are possible. Applications of the approach advocated here together 
with other likelihood methods based on EA are given in [8]. 

Our methodology is a Monte Carlo approach based on two types of proba- 
bilistic constructions. The first of these exploits a duality between diffusions 
and Poisson processes, which has become transparent since the develop- 
ment of EA in [7]. The second involves explicit representations of condi- 
tioned Brownian sample paths constructed so as to be suitably continuous 
in model parameters. Although the entire mathematical construction is quite 
involved, in this paper we are able to distill its implementation down to a 
very simple Monte Carlo algorithm requiring merely collections of Gaussian 
and exponential random variables, as described in Theorem 1. The computer 
code for implementing our algorithm is freely available by e-mail request to 
any of the authors. 

The scope of the methodology we introduce here can be extended in var- 
ious directions. A direct extension is to certain multivariate diffusions. For 
example, it can be seen that our methods extend to processes which (possibly 
after a transformation) can be written as 

dXs=VA{X,;9)ds+ dBs, X^eA'CR™, 

and i? is a m-dimensional Brownian motion. The conditions (Co)-(C3) have 
to be appropriately modified, particularly (C3) now becomes that a = VA 
is such that (a^ + a'){-;9) is bounded above on (zi,oo) x (z2,oo) x ••• x 
{zm, 00) for all 2i, . . . , G R-. Furthermore, it is likely that our methods will 
be extended by currently ongoing work on simulating diffusions for which 
condition (C3) is not required to hold. 

APPENDIX: SIMULATION OF r 

The joint distribution of the minimum and the time instance when this is 
attained (m, r) for a Brownian bridge is of a known form; see [36] . In [7] 
we noticed that the distribution of r conditionally on m can be expressed 
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in terms of the mixture of two inverse Gaussian laws. The exact simulation 
formula is the following: 

= {1 + I[Uo < (1 + • h + I[^7o > (1 + vW^)~'] • l2^}/t, 

where ci = {y — mf' j (2t), C2 = (x — raf' j {2t), and Ii IGau{y^ci/c2, 2ci), 
I2 ~ IGau{\/ C2/C1, 2C2) and Uq ~ Un[0, 1]. We denoted by IGau{-, •) the in- 
verse Gaussian distribution specified by two positive parameters (see, e.g., 
Chapter IV. 4 of [14]). An inverse Gaussian random variable / with param- 
eters {c,d) can be represented as (see page 149 of [14]) 

/ = I[Vo < 1/(1 + z)]cz + I[Vo > 1/(1 + z)]-, 

z 

where z = l + {c/d)Z^/2 - V4(c/d)Z^ + (c/d)^Z4/2, Vq ~ Un[0, 1]. Using this 
formula for Ii, /2 with the same Vq, Z, and after some algebra, the expression 
for r simplifies to (17). 

Acknowledgments. We thank the referees for their valuable suggestions 
and corrections. 
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